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Drawing on similarities in Hartree-Fock theory and the theory of matrix product states (MPS), 
we explore extensions to time evolution, response theory, and Green functions. We derive analytic 
equations of motion for MPS from the least action principle, which describe optimal evolution 
in the small time-step limit. We further show how hnearized equations of motion yield a MPS 
random phase approximation, from which one obtains response functions and excitations. Finally 
we analyze the structure of site-based Green functions associated with MPS, as well as the structure 
of correlations introduced via the fluctuation-dissipation theorem. 



Matrix product states are a powerful class of quantum 
states for one dimensional correlated quantum problems. 
They underpin the density matrix renormalization group 
algorithm which has yielded unprecedented accuracy in 
many systems [ij, 1^ . In an MPS wave function, the quan- 
tum state is approximated by a contraction of indepen- 
dent tensors, each associated with a site on the lattice. 
Because of the contracted product structure, it can be 
viewed as a site-based mean field theory with finite en- 
tanglement. 

A different class of mean field theories is based on in- 
dependent particles. For fermions, this is Hartree-Fock 
(HF) theory, where the quantum state is approximated 
by a Slater determinant of independent electron orbitals. 
HF theory is the foundation on which many other ap- 
proximations are built. For example, time-dependent 
Hartree-Fock theory (TDHF) approximates the true evo- 
lution of a quantum state with a single Slater determi- 
nant. Linearizing the equations of motion yields the ran- 
dom phase approximation (RPA) which gives the exci- 
tation spectrum and response functions l^. Moreover, 
the independent particle picture contained in HF theory 
is the starting point for the Green function approach to 
many particle quantum systems. 

Recently, there has been much work to generalize MPS 
be yon d the ground state formalism to time dependence 
|4l4lO| and response properties llMl4|. Despite the dif- 
ferent physical content of the site-based mean field of the 
MPS and independent particle mean-field of HF, similar 
mathematical structures imply that extensions to MPS 
can be developed in analogy with HF theory. Here, we 
explore this direction. First, we briefly summarize the 
HF and MPS approaches to ground state energies, es- 
tablishing our notation and the parallel structure of the 
theories. We then derive analytic MPS equations of mo- 
tion using the Dirac-Frenkel variational principle. We 
show that they provide an optimal representation of MPS 
time evolution, and discuss their relationship to time evo- 
lution schemes currently in use. Next, we linearize the 
time evolution to derive an MPS analog of the RPA, and 
show how excitation energies and response properties are 
obtained through a linear eigenvalue problem. Finally, 



we explore the structure of Green functions that emerge 
in the MPS theory, and the correlations that arise in MPS 
through the fluctuation-dissipation theorem. 

Stationary States — In HF theory, the A'^-particle wave 
function is approximated by a single Slater determinant 
of orbitals (f>i{r). The variational objects of the theory 
are the orbitals, and the best mean-field approximation 
to the ground state is the set of orbitals that minimizes 
the expectation value of the Hamiltonian, E. This leads 
to the Hartree-Fock equations, f (pi{r) = £i0i(r), where f 
is the Fock operator. Thus, each orbital is an eigenstate 
of a one-particle Hamiltonian that depends on all of the 
orbitals. In this sense, HF theory is a mean field theory 
for independent particles. If the orbitals are expanded 
in a finite basis, then the Hartree-Fock equations may be 
expressed in matrix form: 



F-C==£;SC, 



(1) 



where F is the Fock matrix, 5 is the overlap matrix, and C 
is the coefficient matrix of the orbitals in the finite basis. 
Since F is a function of C, this equation must be solved 
self-consistently. 

The defining equations of an optimal MPS wave func- 
tion can be derived in an analogous way. In an MPS 
wave function for k sites, the amplitude of a Fock state 
is given by the trace of a product of matrices: 



IV')=ETMnAr In), 



(2) 



where |n) = |ni, n2, . . .) labels a state in Fock space. A 
matrix A"' is associated with each site i and orbital occu- 
pancy, Ui. The size of the matrices defines the auxiliary 
dimension M, and the maximum occupation number of 
a site defines the physical dimension d. The d indepen- 
dent M X M matrices associated with each site can be 
collected into a single tensor, (A^)"'^, where a and (3 are 
the auxiliary indices. This tensor can be flattened into a 
vector labeled by a single compound index I — {rii, a, 13). 
Derivatives with respect to the components of the tensor 
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define a (nonorthogonal) basis for the wave function: 

d\i,) 



(3) 



From the linear dependence of \'ip) on each matrix, we 
havelV'HE.Afl^/}.^ 

Requiring the variational energy to be stationary with 
respect to variations in a single tensor defines an eigen- 
value problem: 



Hi ■ — E Si • Aj 



(4) 



The matrices Hi and Si are the Hamiltonian and overlap 
matrices in the local basis defined by A^ . Since the basis 
defined in Eq. (|3]) depends on the other tensors, A; is the 
eigenvector of an effective Hamiltonian defined by the 
other tensors, analogous to Eq. (HI) for the HE orbitals. 

The entire set of equations, one for each tensor, can 
be combined into a single eigenvalue problem by defining 
the compound index fi — {i,ni,a, jS) and collecting all 
the elements of all the tensors into a single vector A. 
This vector contains kdM^ elements. The wave function 
can be expanded as 



(5) 



and the optimal MPS wave function is a self-consistent 
solution of 



H-A = £;S-A. 

The Hamiltonian and overlap matrices are defined as 

1 



h^.--(V'mI^I^-) 



(6) 

(7) 
(8) 



The power of MPS wave functions stems from their 
numerical efficiency. Eor MPS with open boundary con- 
ditions (as used in the DMRG algorithm) expectation 
values of local operators and Hamiltonians can be cal- 
culated with 0{kM^) complexity. The action of local 
operators and Hamiltonians, such as H • A and S • A, 
is also obtained with 0{kM'^) complexity. Solving the 
eigenvalue problem Eq. ([6]) iteratively requires a cost pro- 
portional to the operations H • A and S • A and is thus 
also of 0{kM^) complexity. The pref actor depends on 
preconditioning for both H and S. The DMRG precon- 
ditions S by solving Eq. ^ for a single tensor A^ at a 
time, where the overlap Si can be exactly removed by 
canonicalization. 

Time Evolution — The time-dependent Schrodinger 
equation can be derived by minimizing the Dirac-Frenkel 
action 



dt 



(9) 



with respect to arbitrary variations {5ip\. When the wave 
function is constrained to a particular form, such as a 
Slater determinant or an MPS, variations may only be 
taken with respect to the parameters {Ai} of the wave 
function. 



(10) 



The time derivative of is also constrained in this way. 
Minimizing the action then gives the best approximation 
to the true evolution of the wave function within the space 
of variational wave functions. 

When applied to Hartree-Fock theory, the resulting 
equations of motion are the time-dependent Hartree-Fock 
equations. The time-dependent version of Eq. ([T]) is 



ins-^^ = F{t)-c{t). 



(11) 



The Fock matrix depends on t through the time- 
dependence of the orbitals. It may also depend explicitly 
on t through a time-dependent external potential 

Applying the same variational approach to the action 
of an MPS wave function gives an equation of motion for 
the time dependence of the matrix elements: 



,nS{t)-^^ = H{t)-A{t), 



(12) 



where H and S are the time-dependent versions of Eq. ([7]) 
and Eq. ([Sj. (A similar equation was derived in the 
DMRG context in Ref. See also Ref. Il5i].) Eq. ^ 
and Eq. (fT2|) can be formally solved with a time-ordered 
exponential, which for the MPS takes the form 



Texpj^ / ^"^^^^ ■ ^^^^ 



A(0). (13) 



Moreover, the equations of motion for the MPS can be 
efficiently propagated. For example, if the time interval 
is discretized into units of duration At, then Eq. (|12l) can 
be used to obtain the MPS at tn+i from the MPS at i„. 
Defining e — At/ih, Eq. (fT2|) gives 



S„ • AA„ = e B, 



(14) 



where B„ = H„ ■ A„. This is a linear equation which 
can be solved iteratively with complexity 0{kM^). In 
practice more sophisticated time-propagation schemes, 
such as norm and energy conserving propagators simi- 
lar to those used in time-dependent Hartree-Fock theory, 
should be employed. 

Connections to Time Evolution Algorithms — The 
fundamental difficulty in simulating the evolution of a 
matrix product state is that the matrix dimension re- 
quired to faithfully represent the exact wave function 
grows exponentially in time. The Lagrangian approach 
described above leads to a set of analytic equations for 
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the evolution of an MPS wave function whose auxiliary 
dimension is fixed at M throughout the simulation. In 
contrast, existing methods for the time evolution of MPS 
allow the dimension of the matrix to first grow at each 
time step, then project back onto an MPS of auxihary 
dimension M. Wc now examine the approximate pro- 
jections in existing time evolution algorithms, and show 
that the analytic approach described above is an optimal 
projection in the limit of an infinitesimal time step. 

In time evolution by block decimation (TEBD) 
and the time-dependent density matrix renormalization 
group (t-DMRG) [6-8], the evolution operator exp(eH) 
is factored into a product of local operators using a Trot- 
ter decomposition. Each local evolution operator is then 
applied in sequence. A two-site evolution operator Ui^i+i 
between neighboring sites joins two tensors of the MPS, 
A; and A^+i, into a single object, T^^i+i. The tensor T 
is then approximately projected into a product of two 
tensors of the same size as A^ and A^+i through the sin- 
gular value decomposition. For a single time step (sweep- 
ing over all evolution operators for a local Hamiltonian) 
this update is efficient and is 0{kAI^). However, be- 
cause the approximate projection is only performed lo- 
cally, rather than involving all the tensors, the updated 
matrix product state is not the best representation of 
the evolved wave function for a given auxiliary dimen- 
sion. Furthermore, the algorithm is incompatible with 
long range Hamiltonians. 

In time-dependent matrix product states (tMPS) ,9j, 
the full evolution operator is applied to the current state 
before the projection onto an MPS of auxiliary dimension 
M. In the projection, tMPS attempts to minimize the 
cost function 

A[|V;„+i)] = |||V«+i)-e^«|^„)f (15) 

where \ijjn+i) and \tpn) are the new and old MPS respec- 
tively. The minimization is of complexity 0{kM^), and 
yields in principle the best projected MPS wave func- 
tion. However, since \ipn+i) depends non-linearly on Aj, 
the minimization is not guaranteed to find the optimal 
solution in practice. 

For small time step, however, the projection can be 
done without any non-linear minimization. We see this 
by recognizing that 

^f^^dj^^ , (16) 

ot 

Substituting into the cost function Eq. (|15p and mini- 
mizing with respect to changes in A yields the linear dis- 
cretized equation of motion obtained in Eq. (fT4|) . Thus 
propagation of the MPS equation of motion exactly de- 
termines the optimal projection in a tMPS algorithm, 
without a nonlinear minimization, in the limit At — )■ 0. 

Random Phase Approximation — Excited state prop- 
erties, such as the spectrum and other expectation values. 



can be obtained without studying the full time evolution 
of the system. Only the linearized time evolution, or lin- 
ear response, need be considered. In TDHF, linearization 
of the equations of motion around a stationary state leads 
to the random phase approximation. This is achieved in 
Eq. ([TT|) by taking C{t) = Co + D{t) and expanding all 
quantities to linear order. We now show that a similar 
approach to the MPS equations of motion in Eq. (fT2|) 
yields an MPS analog of the RPA, from which quanti- 
ties such as the excitation spectrum may be obtained. 
These results expand on the analytic response theory we 
recently described in the DMRG context in Ref. [14| . 

We take the zeroth-order MPS to be the ground state, 
defined by A, with energy Eq. The time-dependent MPS 
is defined by 

A{t)=A + h{t). (17) 
The wave function and its derivatives are given by 

m = lJ2[^^+^M\€^) (18) 

I^A.) = lE[A. + b.(i)]IVf.)), (19) 

where |^/;*^°)) and its derivatives are evaluated with b = 0. 
Expanding Eq. (I12p to first order in b gives 

i^iS^^SoS-A + H-b + W-b*. (20) 

The matrices 5 and H are the zeroth-order overlap and 
Hamiltonian matrices. The matrix W couples b to its 
complex conjugate and is defined by 

W^. = i(V'l.°j|H|^(")). (21) 

W is symmetric, but not Hermitian. 

The first term on the right-hand-side of Eq. (pHl) can 
be eliminated by multiplying the entire wave function 1-0) 
by the phase factor ^-^Eot/h^ r^j^^ equations of motion for 
h{t) are harmonic and may be solved by taking h{t) — 
Xe»"t + Y*e~*'^*. Eq. ^ then gives 



"s 







■ H 


w' 







-s* 




w* 


H* 





The eigenvectors of this system of equations define the 
normal modes of the system, and the (positive) eigen- 
values approximate its excitation spectrum. They are 
both efficiently obtained with 0{kM^) complexity. The 
normal modes define the response matrix H(w), which 
determines all the response properties of the system: 
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For example, consider a harmonic perturbation, which 
defines a source q{t) for b: 



(24) 



(The expectation of Q{t) in the ground state is assumed 
to vanish.) Using n(w), the time-dependent variation in 
an observable V due to Q{t) is given by 



■ q(a;) ' 
q*(^) 



(25) 



The linearized equations of motion describe a super- 
position of MPS wave functions. Each tensor product 
Ai ■ A2 . . . An is replaced by a sum of tensor products: 



Ai • A2 . . . Ajv 
Ai • (5A2 . . . Aat + 



5Ai • A2 . . . Aat 

• + Ai ■ A2 . . . SAn 



Thus, the analog of the RPA uses an MPS with auxiliary 
dimension M to give the response properties of an MPS 
of much larger dimension kM. 

Green functions — HP theory is an independent par- 
ticle mean-field theory and the natural Green functions 
that emerge are labeled by particle indices. MPS are 
based on a site-based mean field theory and the natural 
Green functions that emerge are labeled by sites. Here 
we define site-based Green functions, show how they are 
obtained in the RPA formalism above, and analyze the 
nature of the correlations they contain. 

First, we define the one-site density matrix T^^\ which 
is the site-based analog of the one-particle density matrix. 
This is obtained from |^/')('(/'| by tracing out all sites of the 
lattice except site i. The expectation value of a one-site 
operator "P*^'' is given by the trace 



(26) 



Each element Pj*^, is an expectation value of an operator 
%n> defined by 



{n\%n'\n') = SnnSn'n' 



(27) 



For a system with 2 physical degrees of freedom per site, 
the operators a.^ and a| correspond to 701 and 710, and 
the number operator is 711. Creation, annihilation, and 
number operators of systems with more degrees of free- 
dom per site can be constructed from 7 as well. Thus 
the one-site density matrix contains components of one-, 
two-, and mixed-particle density matrices. 

We now define a site-site (retarded) Green function 
n*^*^-'(i). We consider the response of the site density 
matrix at site i and time t, r(*)(i), to a perturbation at 

site J, Q(J) = S{t) Y.nn- ynn'l£,. We then have 



n. 



(y) 



nn :nn' 



(28) 



Note that by choosing appropriate combinations of 7,*^, 
and 7^>j/, it is possible to construct spectral functions 
and one-, two-particle, and other types of conventional 
Green functions. 

Within the RPA framework, the elements of II^'-'^ are 
obtained from the MPS response matrix 11^,^ as 

^',nn'^) = (^l7i^-IV^.). (29) 

The response function and the correlation functions 
of the ground state are related through the fiuctuation- 
dissipation theorem 3, 12 [• This allows us to analyze 
the nature of the additional correlations introduced into 
the MPS at the RPA level. We have 



\ Inn' Inn' 



1 

2^ 



cf^n(;^)_,(w) (30) 



From the expression for the response matrix Ii^i,{uj) 
above, we see that the correlation function is composed 
of a sum of terms, each of the form 

= (V^|7i^,-Axx-7itlV') (31) 



plus the corresponding contributions from Kxy ^ Ayx, 
and Ayy from Eq. (j23p . Axx is a matrix product oper- 
ator (MPO), and is responsible for additional correlation 
introduced at the RPA level relative to the ground state 
(which is obtained by setting all A = 1). By regrouping 
terms in the sum, we obtain 



(32) 



Each term in parentheses is a sum of k MPS of dimen- 
sion M, and their outer product is formally an MPO of 
dimension 2kM. The sum includes an MPO for each 
normal mode, making Axx an MPO of very large auxil- 
iary dimension. Thus, the entanglement and correlation 
possible at the RPA level are greater than the auxiliary 
dimension M of the matrices might suggest. 

Conclusions — In this paper, we have drawn on par- 
allels between the product structure of HF theory and 
MPS to develop analytic equations of motion for time 
evolution, an MPS RPA for response and excitations, and 
a theory of site-based Greens functions which introduce 
correlations beyond the MPS ground state, as demon- 
strated via the fiuctuation-dissipation theorem. The new 
time evolution and response algorithms proposed here 
can all be implemented with a complexity proportional 
to the ground state DMRG algorithm. Finally, our work 
suggests that further extensions of Hartree-Fock theory 
may also be carried over to MPS, and these will be ex- 
plored in the future. 



5 



This work was supported by the the Cornell Center 
for Materials Research, the Center for Molecular Interfac- 
ing, NSF CAREER, DOE CSGF, the Camille and Henry 
Dreyfus Foundation, the David and Lucile Packard Foun- 
dation, and the Alfred P. Sloan Foundation. 



* Electronic address: 'gc238@cornell.edu' 
[1] S. White, Phys. Rev. B 48, 10345 (1993). 
[2] S. Rommer and S. Ostlund, Phys. Rev. B 55, 2164 
(1997). 

[3] A. McLachlan and M. Ball, Rev. Mod. Phys. 36, 844 
(1964). 

[4] H. G. Luo, T. Xiang, and X. Q. Wang, Phys. Rev. Lett. 

91, 049701 (2003). 
[5] G. Vidal, Phys. Rev. Lett. 93, 40502 (2004). 



[6] S. White and A. Feiguin, Phys. Rev. Lett. 93, 76401 
(2004). 

[7] A. Daley C. KoUath, U. SchoUwock, and G. Vidal, J. 

Stat. Mech. 2004, P04005 (2004). 
[8] P. Schmitteckert, Phys. Rev. B 70, 121302(R) (2004). 
[9] F. Verstraete, J. Garcia-RipoU, and J. Cirac, Phys. Rev. 

Lett. 93, 207204 (2004). 
[10] U. SchoUwock, Ann. Phys. 326, 96 (2010). 
[11] K. Hallberg, Phys. Rev. B 52, 9827 (1995). 
[12] T. Kiihner and S. White, Phys. Rev. B 60, 335 (1999). 
[13] E. Jeckelmann, Phys. Rev. B 66, 45114 (2002). 
[14] J. Dorando, J. Hachmann, and G. Chan, J. Chem. Phys. 

130, 184111 (2009). 
[15] K. Ueda, C. Jin, N. Shibata, Y. Hieida, and T. Nishino, 

arXiv:cond-mat/0612480, (2006). 
[16] H. Callen and T. Welton, Phys. Rev. 83, 34 (1951). 
[17] D. Pines and P. Nozieres, Theory of Quantum Liquids: 

Normal Fermi Liquids (Westview Press, 1994). 



